#include "mt3d.h"

void mt3d::save_model(std::string filename, bool debug)
{
    gctl::gmshio modelio;

    std::clog << "Writing Gmsh file ... ";
    modelio.init_file(filename, gctl::Output);
    modelio.set_packed(gctl::NotPacked, gctl::Output);
    modelio.save_mesh(mesh_elems_, mesh_nodes_);
    modelio.save_data("Node tag", node_markers_, gctl::NodeData);
    modelio.save_data("Element tag", elem_markers_, gctl::ElemData);

    if (debug) // save electronic field
    {
        gctl::array<double> node_vol(node_num_, 0.0);
        gctl::array<double> te_xval_real(node_num_, 0.0), te_xval_imag(node_num_, 0.0);
        gctl::array<double> te_yval_real(node_num_, 0.0), te_yval_imag(node_num_, 0.0);
        gctl::array<double> te_zval_real(node_num_, 0.0), te_zval_imag(node_num_, 0.0);
        gctl::array<double> tm_xval_real(node_num_, 0.0), tm_xval_imag(node_num_, 0.0);
        gctl::array<double> tm_yval_real(node_num_, 0.0), tm_yval_imag(node_num_, 0.0);
        gctl::array<double> tm_zval_real(node_num_, 0.0), tm_zval_imag(node_num_, 0.0);

        gctl::vertex3dc* tmp_vert[4];
        double x, y, z, ux, uy, uz, dummy;
        complex_d te_fx, te_fy, te_fz, tm_fx, tm_fy, tm_fz;
        efem_sf::ordertype_e o_type;

        double tmp_vol, elem_cndt;
        for (size_t e = 0; e < elem_num_; e++)
        {
            // 双重加权平均 同时考虑三角形面积与电导率的权重 减少空气层界面的畸变
            get_element_conductivity(elem_markers_[e], elem_cndt);
            tmp_vol = mesh_elems_[e].volume() * elem_cndt;

            for (size_t n = 0; n < 4; n++)
            {
                node_vol[mesh_elems_[e].vert[n]->id] += tmp_vol;
            }
        }

        for (size_t e = 0; e < elem_num_; e++)
        {
            tmp_vert[0] = mesh_elems_[e].vert[0];
            tmp_vert[1] = mesh_elems_[e].vert[1];
            tmp_vert[2] = mesh_elems_[e].vert[2];
            tmp_vert[3] = mesh_elems_[e].vert[3];

            get_element_conductivity(elem_markers_[e], elem_cndt);
            tmp_vol = mesh_elems_[e].volume() * elem_cndt;

            for (size_t n = 0; n < 4; n++)
            {
                x = mesh_elems_[e].vert[n]->x;
                y = mesh_elems_[e].vert[n]->y;
                z = mesh_elems_[e].vert[n]->z;

                te_fx = te_fy = te_fz = tm_fx = tm_fy = tm_fz = _zero;
                for (size_t i = 0; i < 6; i++)
                {
                    if (elem_edge_reversed_[e][i]) o_type = efem_sf::Reverse;
                    else o_type = efem_sf::Normal;

                    elsf_.tetrahedron(x, y, z, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, tmp_vert[3]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, tmp_vert[3]->y, 
                        tmp_vert[0]->z, tmp_vert[1]->z, tmp_vert[2]->z, tmp_vert[3]->z, 
                        i, efem_sf::Value, o_type, ux, uy, uz, dummy, dummy, dummy);

                    te_fx += ux*X_TE_[elem_edge_index_[e][i]];
                    te_fy += uy*X_TE_[elem_edge_index_[e][i]];
                    te_fz += uz*X_TE_[elem_edge_index_[e][i]];
                    tm_fx += ux*X_TM_[elem_edge_index_[e][i]];
                    tm_fy += uy*X_TM_[elem_edge_index_[e][i]];
                    tm_fz += uz*X_TM_[elem_edge_index_[e][i]];
                }

                te_xval_real[mesh_elems_[e].vert[n]->id] += te_fx.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_xval_imag[mesh_elems_[e].vert[n]->id] += te_fx.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_yval_real[mesh_elems_[e].vert[n]->id] += te_fy.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_yval_imag[mesh_elems_[e].vert[n]->id] += te_fy.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_zval_real[mesh_elems_[e].vert[n]->id] += te_fz.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_zval_imag[mesh_elems_[e].vert[n]->id] += te_fz.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];

                tm_xval_real[mesh_elems_[e].vert[n]->id] += tm_fx.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_xval_imag[mesh_elems_[e].vert[n]->id] += tm_fx.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_yval_real[mesh_elems_[e].vert[n]->id] += tm_fy.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_yval_imag[mesh_elems_[e].vert[n]->id] += tm_fy.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_zval_real[mesh_elems_[e].vert[n]->id] += tm_fz.real()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_zval_imag[mesh_elems_[e].vert[n]->id] += tm_fz.imag()*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
            }
        }

        modelio.save_data("TE_Ex_real", te_xval_real, gctl::NodeData);
        modelio.save_data("TE_Ex_imag", te_xval_imag, gctl::NodeData);
        modelio.save_data("TE_Ey_real", te_yval_real, gctl::NodeData);
        modelio.save_data("TE_Ey_imag", te_yval_imag, gctl::NodeData);
        modelio.save_data("TE_Ez_real", te_zval_real, gctl::NodeData);
        modelio.save_data("TE_Ez_imag", te_zval_imag, gctl::NodeData);
        modelio.save_data("TM_Ex_real", tm_xval_real, gctl::NodeData);
        modelio.save_data("TM_Ex_imag", tm_xval_imag, gctl::NodeData);
        modelio.save_data("TM_Ey_real", tm_yval_real, gctl::NodeData);
        modelio.save_data("TM_Ey_imag", tm_yval_imag, gctl::NodeData);
        modelio.save_data("TM_Ez_real", tm_zval_real, gctl::NodeData);
        modelio.save_data("TM_Ez_imag", tm_zval_imag, gctl::NodeData);

        /*
        gctl::array<double> te_xval_real(elem_num_, 0.0), te_xval_imag(elem_num_, 0.0);
        gctl::array<double> te_yval_real(elem_num_, 0.0), te_yval_imag(elem_num_, 0.0);
        gctl::array<double> te_zval_real(elem_num_, 0.0), te_zval_imag(elem_num_, 0.0);
        gctl::array<double> tm_xval_real(elem_num_, 0.0), tm_xval_imag(elem_num_, 0.0);
        gctl::array<double> tm_yval_real(elem_num_, 0.0), tm_yval_imag(elem_num_, 0.0);
        gctl::array<double> tm_zval_real(elem_num_, 0.0), tm_zval_imag(elem_num_, 0.0);

        gctl::vertex3dc* tmp_vert[4];
        double x, y, z, ux, uy, uz, dummy;
        complex_d te_fx, te_fy, te_fz, tm_fx, tm_fy, tm_fz;
        efem_sf::ordertype_e o_type;

        for (size_t e = 0; e < elem_num_; e++)
        {
            tmp_vert[0] = mesh_elems_[e].vert[0];
            tmp_vert[1] = mesh_elems_[e].vert[1];
            tmp_vert[2] = mesh_elems_[e].vert[2];
            tmp_vert[3] = mesh_elems_[e].vert[3];

            x = 0.25*(tmp_vert[0]->x + tmp_vert[1]->x + tmp_vert[2]->x + tmp_vert[3]->x);
            y = 0.25*(tmp_vert[0]->y + tmp_vert[1]->y + tmp_vert[2]->y + tmp_vert[3]->y);
            z = 0.25*(tmp_vert[0]->z + tmp_vert[1]->z + tmp_vert[2]->z + tmp_vert[3]->z);

            te_fx = te_fy = te_fz = tm_fx = tm_fy = tm_fz = _zero;
            for (size_t i = 0; i < 6; i++)
            {
                if (elem_edge_reversed_[e][i]) o_type = efem_sf::Reverse;
                else o_type = efem_sf::Normal;

                elsf_.tetrahedron(x, y, z, 
                    tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, tmp_vert[3]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, tmp_vert[3]->y, 
                    tmp_vert[0]->z, tmp_vert[1]->z, tmp_vert[2]->z, tmp_vert[3]->z, 
                    i, efem_sf::Value, o_type, ux, uy, uz, dummy, dummy, dummy);

                te_fx += ux*X_TE_[elem_edge_index_[e][i]];
                te_fy += uy*X_TE_[elem_edge_index_[e][i]];
                te_fz += uz*X_TE_[elem_edge_index_[e][i]];
                tm_fx += ux*X_TM_[elem_edge_index_[e][i]];
                tm_fy += uy*X_TM_[elem_edge_index_[e][i]];
                tm_fz += uz*X_TM_[elem_edge_index_[e][i]];
            }

            te_xval_real[e] = te_fx.real();
            te_xval_imag[e] = te_fx.imag();
            te_yval_real[e] = te_fy.real();
            te_yval_imag[e] = te_fy.imag();
            te_zval_real[e] = te_fz.real();
            te_zval_imag[e] = te_fz.imag();

            tm_xval_real[e] = tm_fx.real();
            tm_xval_imag[e] = tm_fx.imag();
            tm_yval_real[e] = tm_fy.real();
            tm_yval_imag[e] = tm_fy.imag();
            tm_zval_real[e] = tm_fz.real();
            tm_zval_imag[e] = tm_fz.imag();
        }

        modelio.save_data("TE_Ex_real", te_xval_real, gctl::ElemData);
        modelio.save_data("TE_Ex_imag", te_xval_imag, gctl::ElemData);
        modelio.save_data("TE_Ey_real", te_yval_real, gctl::ElemData);
        modelio.save_data("TE_Ey_imag", te_yval_imag, gctl::ElemData);
        modelio.save_data("TE_Ez_real", te_zval_real, gctl::ElemData);
        modelio.save_data("TE_Ez_imag", te_zval_imag, gctl::ElemData);
        modelio.save_data("TM_Ex_real", tm_xval_real, gctl::ElemData);
        modelio.save_data("TM_Ex_imag", tm_xval_imag, gctl::ElemData);
        modelio.save_data("TM_Ey_real", tm_yval_real, gctl::ElemData);
        modelio.save_data("TM_Ey_imag", tm_yval_imag, gctl::ElemData);
        modelio.save_data("TM_Ez_real", tm_zval_real, gctl::ElemData);
        modelio.save_data("TM_Ez_imag", tm_zval_imag, gctl::ElemData);
        */
    }

    std::clog << "done\n";
    return;
}
